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This research combines Volterra theory and proper orthogonal decomposition (POD) 
into a hybrid methodology for reduced-order modeling of aeroelastic systems. The out- 
come of the method is a set of linear ordinary differential equations (ODEs) describing 
the modal amplitudes associated with both the structural modes and the POD basis 
functions for the fluid. For this research, the structural modes are sine waves of varying 
frequency, and the Volterra-POD approach is applied to the fluid dynamics equations. 
The structural modes are treated as forcing terms which are impulsed as part of the 
fluid model realization. Using this approach, structural and fluid operators are coupled 
into a single aeroelastic operator. This coupling converts a free boundary fluid problem 
into an initial value problem, while preserving the parameter (or parameters) of interest 
for sensitivity analysis. The approach is applied to an elastic panel in supersonic cross 
flow. The hybrid Volterra-POD approach provides a low-order fluid model in state-space 
form. The linear fluid model is tightly coupled with a nonlinear panel model using an 
implicit integration scheme. The resulting aeroelastic model provides correct limit-cycle 
oscillation prediction over a wide range of panel dynamic pressure values. Time inte- 
gration of the reduced-order aeroelastic model is four orders of magnitude faster than 
the high-order solution procedure developed for this research using traditional fluid and 
structural solvers. 


Introduction 

Volterra methods 1 * and proper orthogonal decompo- 
sition 2,3 (POD) are two of the more prevalent reduced- 
order modeling (ROM) techniques well-suited to non- 
linear dynamics. 4-8 The application of ROM tech- 
niques to aeroelastic systems is an active area of re- 
search, motivated by the desire for faster algorithms 
that are well-suited to the design environment for air- 
craft. For example, transonic, fluid-structure inter- 
action is a particular application of interest to both 
external and internal aerodynamicists because moving 
shock waves in the flow necessitate high-fidelity numer- 
ical flow solvers which are too cumbersome for iterative 
design analysis. Regardless of the application, when 
nonlinearities are present in either the flow field or the 
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structure, established order-reduction methods that 
rely on linearized dynamics are of little use. 

Over the past three years, applications of POD 
to the Euler equations have produced reduced order 
aeroelastic models that properly capture aerodynamic 
nonlinearities. A low-order POD representation of 
the discrete, 2-D Euler equations 9 was coupled with 
the von Karman equation to simulate the dynamics 
of flow over a flexible panel. 10 Subsequently, a new 
approach was taken, involving domain decomposition, 
that allowed LCO to be accurately simulated in the 
transonic regime. 11 In that study, full-order and re- 
duced order models of a small flow region containing 
a moving shock were decomposed from the larger flow 
domain. Both approaches enabled a physically consis- 
tent treatment of the aerodynamic nonlinearity. In a 
more recent paper, 12 the original POD/ROM method- 
ology used for flow over an elastic panel 10 was revis- 
ited to improve the temporal coupling between the 
aerodynamic and structural dynamic equations. Fur- 
thermore, a modal representation of the structure was 
employed, which permitted a more efficient formula- 
tion of the reduced-order aeroelastic system. 

All of the studies mentioned above relied on a ROM 
technique called subspace projection for time integra- 
tion of the reduced-order model. While sufficient to 
demonstrate the accuracy of the POD basis functions, 


1 of 12 


American Institute of Aeronautics and Astronautics Paper 2003-1922 



subspace projection was not an efficient way to time 
integrate the low-order, aeroelastic ROM. Generally, 
four orders of magnitude reduction in fluid system de- 
grees of freedom (DOFs) were demonstrated in the 
above studies. Time integrating these POD/ROMs 
with subspace projection generally produced about 
one order of magnitude improvement in compute time 
to accompany a much larger drop in problem order. 

The applicability of POD basis functions to nonlin- 
ear problems has been documented in the literature, 
but a tractable nonlinear, low-order model realiza- 
tion procedure is a key missing link. Two techniques, 
Galerkin projection and direct projection, have been 
recently reported as having potential for obtaining 
nonlinear terms for POD/ROMs. 13 However, the lin- 
ear portion of these realization procedures is generally 
unstable, requiring dissipation techniques that affect 
model performance. The Volterra-POD approach pro- 
vides a stable reduced-order equation set, and is an 
important advance toward achieving stable, nonlinear 
reduced-order models. 

The hybrid Volterra-POD method was recently de- 
veloped to replace subspace projection for time inte- 
gration of POD/ROMs applied to compressible flow 
fields. 14 The goal of the Volterra-POD approach was 
to achieve computational savings on the order of DOF 
reductions. This goal was achieved in the initial ap- 
plication, where four orders of magnitude reduction 
was obtained in both DOFs and compute time. To 
date, the hybrid Volterra-POD method has only been 
applied to subsonic flow-fields characterized by linear 
behavior, with fixed boundaries. The product of the 
technique was a linear, state-space system of ODEs 
governing the dynamics of modal coefficients corre- 
sponding to a small number of POD basis functions. 
The state-space realization was obtained from a set of 
impulse responses that were processed using the Eigen- 
system Realization Algorithm (ERA). 15 ’ 16 

This research will extend the Volterra-POD ap- 
proach to supersonic flow-fields with dynamic bound- 
ary behavior. The POD-Volterra method will be ap- 
plied to a two-dimensional elastic panel in inviscid, 
supersonic cross-flow. The Volterra-POD approach 
will be used to identify a low-dimensional, linear 
POD/ROM for the fluid. The POD/ROM will be 
tightly coupled to a low-dimensional, nonlinear model 
of the von Karman plate equation. 17 The aeroelas- 
tic response will be obtained using an implicit time- 
integration scheme. 

The Volterra-POD technique involves procedures 
that require the selection of parameters such as im- 
pulse size, several data windowing lengths, and im- 
pulse sampling frequency. The choice of POD basis 
affects performance as well. Some considerations for 
generating the POD basis include choice of base flow 
(the POD/ROM determines the perturbation to this 
base flow), snapshot collection method and sampling 


frequency (the method of snapshots 18 will be discussed 
in the next section). 

The research will consider two base flow cases, and 
two snapshot collection methods. Both uniform flow 
at free stream conditions, and steady-state flow over 
a static panel deflection will be considered as base 
flow cases. An aeroelastic POD basis will be gener- 
ated by sampling a small portion of the time history 
for a baseline LCO case, which was the approach in 
recent applications using subspace projection for this 
problem. 12 In addition, we will investigate using the 
impulse response of the fluid system to generate a POD 
basis. The full-system impulse response is collected as 
part of the Volterra-POD approach, and the impulse 
responses can be sampled as snapshots in lieu of the 
LCO time-history. Finally, we will apply POD to the 
structural dynamics, couple the structural POD/ROM 
with the fluid POD/ROM and examine performance. 
We will record the various parameter settings used to 
generate the aeroelastic POD/ROMs for each case. 

The linearity of the supersonic flow-field will be ex- 
amined as part of the ROM analysis. The principle 
of superposition applies in a linear flow- field, which 
enables a host of linear order-reduction techniques, in- 
cluding the Volterra-POD technique detailed in this 
paper. While the supersonic, aeroelastic flow-field is 
well represented by a linear fluid model, we will demon- 
strate that the supersonic flow-field itself is not linear 
in general. 

The performance of the Volterra-POD aeroelastic 
ROMs will be quantified in accuracy, order reduction, 
and computational savings. A high-order, full-system 
representation of the problem is required for snapshot 
collection. The flow held and panel response for the 
full-system model will serve as the baseline for per- 
formance comparison. Accuracy will be quantified by 
comparing LCO panel response, how-held pressure dis- 
tribution on the elastic panel, LCO frequency, and 
LCO phase for a variety of panel dynamic pressure 
values. Finally, the robustness of the Volterra-POD 
method for predicting LCO response across a broad 
parameter space will also be addressed. 

Formulation 

This section describes the full-system aeroelastic 
model, introduces POD, and overviews Volterra meth- 
ods. In addition, we fully develop the Volterra-POD 
approach and the synthesis of aeroelastic ROMs. 

Structural Dynamics Equations 

Two-dimensional how over a semi-infinite, pinned 
panel of length L is considered. Panel dynamics are 
computed with von Karman’s large-dehection plate 
equation, which is placed in nondimensional form us- 
ing aerodynamic scales L and Woo 19 (0 < x < 1): 
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dw 


(l-, 2 )/ dx. (2) 


The nonlinear, in-plane load in Eqn. (2), serves to limit 
panel deflections w(x,t ) induced by fluid-structure in- 
teraction. Here, the load is assumed to be distributed 
uniformally over the panel. 17 Equation (1) is compa- 
rable to similar formulations found in the literature, 17 
although Wd and t,j are scaled by h and (p s ft,L 4 ) , 
respectively. Two pinned boundary conditions are en- 

o2 

forced at the panel’s endpoints: w = 0 and = 0. 

A modal solution for the deflection w( x, t) is as- 
sumed: 

m s 

w(x, t) = '^2 cLi{t) sin(i7ra;), (3) 

i= 1 

where m s is the number of structural modes retained, 
and the modal amplitudes a.; vary in time and are col- 
located in the array a. The Galerkin method is used 
to obtain a low-order set of ordinary-differential equa- 
tions describing the behavior of aj. 17 First, Eqn. (3) 
is substituted into Eqn. (1). The resulting expres- 
sion is then integrated, following pre-multiplication by 
sin(wrx), to yield ( i = 1, ...,m s ) 


1 .. /i(i7r) 4 6/u / h 


2 a,+ 2A 


-T il p-o 


2\ ( 77r ) 2 D 

1 a—z—di = uPf 


where a = J2 r a: 


2 (rir) 2 


( 4 ) 


P, = 


r r 2 

1 ( i 


and 


V7 Ml 


— p sin(inx)dx. (5) 


The projected pressure components, Pi, are integrated 
from the aerodynamic solution with the midpoint rule, 
using flowfield pressures obtained at grid points on the 
panel surface. 12 The aerodynamic equations, their dis- 
cretization, and their solution are discussed in later 
sections. While equivalent to other formulations in 
the literature, Eq. (4) has two notable distinctions. 
First, the different form of scaling described above al- 
ters equation coefficients, and, second, an expression 
relating p to the state of the panel is not assumed. 12 

The structural dynamics equation Eqn. (4) is placed 
in first-order form by introducing a mode speed array, 
b, such that a, = 6,, 


Fluid Dynamics 

The dynamics of inviscid fluid flows are governed by 
the Euler equations. The two-dimensional Euler equa- 
tions are given below in strong conservation form: 20 


dU 

~dt 


U(x, t) 


dE(U) dF(U) 
dx dy 

P 

m x 

m. y ’ 

E t 


(8a) 

(8b) 


where p, m x , m y and Et are functions of space and 
time. Since we assume an ideal gas for our applica- 
tions, this equation set can be closed using the ideal 
gas law. 

The solution of the Euler equations can be approx- 
imated using either finite-difference, finite- volume, or 
finite-element techniques. To do this, the spatial do- 
main is discretized, and the flow variables in U(x,t) at 
each discrete location are collocated into a column vec- 
tor U(t). Time integration across the computational 
mesh is used to obtain flow solutions. 

Since the Euler equations are linear in the time 
derivative, and quasi-linear in the spatial deriva- 
tive, 20,21 the spatial derivatives and the time deriva- 
tives in Eqn. (8a) can be separated to form an evolu- 
tionary system. To accomplish this, the spatial deriva- 
tives of the flux terms ^ and ^ are grouped to form 
a nonlinear operator R acting on the set of fluid vari- 
ables. The fluid dynamics from Eqn. (8a) can then be 
expressed as 


dU(x,t) 

dt 


R(U(x,t)) . 


(9) 


When discretized this expression takes the form 


^ = R(U(t)) ■ ( 10 ) 


Equation (10) is referred to as the full-system dynam- 
ics. 

A finite-volume scheme was the basis for the full- 
order solver used in this research, which approximated 
the integral form of the Euler equations: 

f UdV + f ( Ei+ Fj) ■ dS = 0 . (11) 

dt Jy J Qy 


k = 


p(^) 4 

X 


6/i / Lin 

T l ~h 


(1 - , 2 ) , 


; + 2 pPi- 


( 6 ) 

The mode speeds and amplitudes are collocated into 
a structural solution array, Y s , leading to a general 
form of the structural equation: 


Y s = [b, a] T (7a) 

Y s = R S (Y S ,P- p, A, h/L). (7b) 


The grid points in the computational mesh described 
earlier were used to form corners for cells. For each 
cell, the integral form of the Euler equations reduced 
to the following, assuming no grid deformation: 


d_ 

dt 


Ui,j + ^ + Pi,jj) ' 

sides 


dSij 

dAi,j 


= 0 . 


(12) 


The flux terms Eij and F- t j were computed using 
second-order Roe averaging, 20 and the flow variables 
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Uij were evaluated as cell averages. Time integration 
across the computational mesh was used to obtain flow 
solutions. This was accomplished with a first-order- 
accurate, forward Euler approximation. 

External boundaries were handled with ghost cells. 
The fluid values for the ghost cells at the far 
field boundaries were determined using characteristic 
boundary conditions. 20 The bump-surface was mod- 
eled using a transpiration approximation. 22 The finite- 
volume fluid solver and the transpiration boundary 
condition were validated using a combination of the- 
ory and experimental data. Subsonic performance was 
validated using wind-tunnel data. 23 Supersonic per- 
formance was validated using oblique shock theory. 24 
Time-accurate performance was validated by correctly 
predicting the vortex shedding frequency from a cylin- 
der in cross flow. 

Time Integration of the Coupled Full-Order 
Equations 

The systems of discretized fluid dynamic equations, 
U(t), and modal structural equations, Y a , are com- 
bined into a single time-dependent system represen- 
tative of the complete interaction between structure 
and inviscid flow. Time integration proceeds in two 
steps, assuming an O(At) lag in the synchronization 
of fluid and structure. First, the structural variables 
are updated from time level n to n + 1 using a Crank- 
Nicolson procedure to be described below (but limited 
here to only structural variables). During this step, the 
pressures known at grid points on the panel surface are 
considered frozen. In the second step, the aerodynamic 
variables are explicitly updated using only structural 
variables defined at time level n. 

Grid Generation and Time Step 

The flow is simulated over a physical domain of 
length Dl, centered about x = 0, and height Dh ■ The 
domain is discretized using I nodes in the streamwise 
direction and J nodes normal to the panel. Indices i 
(1 < i < I) and j (1 < j < J) are used to denote grid 
points comprising corners of cells for the finite-volume 
scheme. Grid points are clustered in the direction nor- 
mal to the panel at the panel surface, with minimum 
spacing denoted by A wa u. The spacing of grid points is 
specified to grow geometrically with j from the panel 
boundary. In the streamwise direction, the node spac- 
ing is chosen to be uniform over the deforming panel 
segment (coincident with the structural grid), while 
growing geometrically upstream of the leading edge 
(positioned at * = Ile) and downstream of the trail- 
ing edge (positioned at i = Ite )• Calculations are 
carried out with a baseline grid given by the following: 
I = 141, J = 116, D l = 50, D h = 25, I LE = 45, 
Ite — 97, and A wa n = 0.0125. 


Proper Orthogonal Decomposition 

POD is a technique to identify a small number of 
basis functions that adequately describe the behav- 
ior of the full-system dynamics (Eqn. (10)) across 
some parameter space of interest. A summary of POD 
as it applies to a spatially-discretized flow field fol- 
lows. A detailed description of POD is available in the 
literature. 4,8 For simplicity, consider only one fluid 
variable, w which when spatially descritized us- 
ing N nodes is denoted w(f). For this fluid variable, 
the full-system dynamics in Eqn. (10) is expressed as 

^ = R w ( w) . (13) 

Spectral methods approximate the solution w (X,t) as 

M 

w(®, t) « ^2 a k (t)(j) k (x ) . (14) 

k = 1 

When the domain is spatially discretized, </>k(x) be- 
comes a vector cf> k , and the following relation applies: 

M 

( 15 ) 

fc = 1 

The set of vectors {<j> fc } are discrete basis functions 
corresponding to the computational mesh defined for 
the numerical solver. The set {a^} are the modal coef- 
ficients, and Eqn. (15) can be represented using matrix 
algebra. The fluid modes comprise columns of a modal 
matrix 'I>, and the coefficients are collocated into a 
column vector w (f). POD produces a linear transfor- 
mation $ between the full-order solution, w, and the 
reduced-order solution, w: 

w (t) = W 0 + 4>w(t) . (16) 

The reduced-order variable w(t) represents deviations 
of w (t) from a base solution Wo- The subtraction of 
Wo will result in zero- valued boundaries for the POD 
modes wherever constant boundary conditions occur 
on the domain. 

$ is constructed by collecting observations of the 
solution w(1)-Wq at different time intervals through- 
out the time integration of the full-system dynamics. 
These observations are called snapshots 18 and are gen- 
erally collected to provide a good variety of flow field 
dynamics while minimizing linear dependence. The 
snapshot generation procedure is sometimes referred 
to as POD training. 8 

A total of Q snapshots are collected from the full- 
system dynamics. These are vectors of length N. The 
set of snapshots describe a linear space that is used 
to approximate both the domain and the range of the 
nonlinear operator R w . The linear space is defined 
by the span of the snapshots. 25 POD identifies a new 
basis for this linear space that is optimally convergent 4 
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in the sense that no other set of basis functions will 
capture as much energy in as few dimensions as the 
POD basis functions. To identify the POD basis, the 
snapshots are compiled into an N x Q matrix S', known 
as the snapshot matrix. The mapping function d> is 
then developed using 

S T SV = VA , (17a) 

= SV . (17b) 

Here V is the matrix of eigenvectors of S T S , and A is 
the corresponding diagonal matrix of eigenvalues. To 
eliminate redundancy in the snapshots, the columns 
of V corresponding to very small eigenvalues in A are 
truncated. The matrix of eigenvalues A is also resized 
to eliminate the rows and columns corresponding to 
the removed eigenvalues. If Q — M columns of V are 
truncated, the resulting reduced order mapping <I> will 
be an N x M matrix. The reduced-order mapping 
represents an isomorphism between N and M dimen- 
sional space. It determines the coordinates of w (t) in 
terms of the M remaining basis functions, cf> k . 

The reduced-order mappings for each fluid variable 
are developed separately, and individual S and V ar- 
rays are collocated as blocks into a larger set of arrays, 
also denoted S and V, to form 

U(t) « U 0 + $U(t) , (18a) 

$ = SV . (18b) 

These versions of Eqn.’s (16) and (17b), respectively, 
apply to the entire set of fluid variables. 

POD of the discrete, panel position vector w(x, t ) — > 
w(t) and panel velocity vector s(t) = w(t) is accom- 
plished in a similar manner as described above for the 
fluid system. Unlike the fluid POD basis functions, 
there is no base term subtracted from the snapshots 
when generating a structural POD basis. 

Volterra Methods 

Consider time- invariant, nonlinear, continuous- 
time, systems. Of interest is the response of the system 


about an initial state X(0) = Xq due to an arbi- 
trary input u(t) (we take u as a real, scalar input) 
for t > 0. As applied to these systems, Volterra the- 
ory 1, 26-28 yields the response 

X(t) 

= h 0 + f h\(t — t)u(t)cIt (19) 

Jo 

+ 

/ / h 2 (t - Ti,t - T 2 )u{T 1 )u(T 2 )dT 1 dT 2 

Jo Jo 

+ 

N ,.t ft. 

^ ^ / " 1 ^ n (J l”l , ■ t T n ) 

n —3 J ° J 0 


u(n).. u(T n )dTi.. dr n . 


formulation is considered: 

X(t) = ho+ f h\{t — t)u{t)<It 
Jo 

+ [ f h 2 (t — Ti,t — T 2 )u(Ti)u(T2)dTidT2 . (20) 

Jo Jo 

The assumption of a weakly nonlinear system is con- 
sistent with the emergence of limit-cycle oscillation of 
a 2-D aeroelastic system in transonic flow through a 
supercritical Hopf bifurcation. 29 For linear systems, 
only the first-order kernel is non-trivial, and there are 
no limitations on input amplitude. 

The first- and second-order kernels are presented be- 
low in final form: 5 

Mn) = 2X 0 (t 1 )-- 2 X 2 (t 1 ) , (21) 

Mn,r 2 ) = ^(X 1 (r 1 ,r 2 )-X o (r 1 )-X o (r 2 )022) 

In (21), X 0 (ti) is the time response of the system to 
a unit impulse applied at time 0 and X 2 (ti) is the 
time response of the system to an impulse of twice 
unit magnitude at time 0. These response functions 
represent the memory of the system. If the system is 
linear, then X 2 = 2Xq and h i = Xq, which is why 
the first-order kernel is referred to as the linear unit 
impulse response. The identification of the second- 
order kernel is more demanding, since it is dependent 
on two parameters. Assuming r 2 > t\ in (22), X 0 (t 2 ) 
is the response of the system to an impulse at time r 2 . 

Time is discretized with a set of time steps of equiv- 
alent size. Time levels are indexed from 0 (time 0) to 
n (time £), and the evaluation of X at time level n is 
denoted by X[n]. The convolution in discrete time is 

N 

X[n] = ho + h\]n — fc]w[fc] (23) 

fc= o 
N N 

+ ^2 h 2 [n — k\, n — k 2 \u[k\\u[k 2 ](2A) 

ki—0 /c2— 0 

The linearized and nonlinear Volterra kernels can 
be transformed into linearized and nonlinear (bilinear) 
state-space systems that can be easily implemented 
into other disciplines such as controls and optimiza- 
tion. 5,27 For linear dynamics, state-space realization 
using the Eigensystem Realization Algorithm (ERA) 
has been used to generate linear, aeroelastic systems. 16 
Nonlinear system realization is an active area of re- 
search. 

System Realization 

The ERA method 15 identifies a discrete, linear, 
time-invariant state-space realization of the form, 


The Volterra series can be accurately truncated be- 
yond the second-order term when a weakly nonlinear 


X[n +1] = AX[n\ + Bu[n\ 

Y[n] = CX[n\ , (25) 
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using data from a complete ensemble of impulse re- 
sponses. Initial state responses can be used in lieu 
of impulse responses, but we only consider impulse 
response data in this overview for simplicity. The sys- 
tems realization procedure takes measurement data 
Y[n] from the free response of the system and pro- 
duces a minimal state-space model A , B , and C such 
that functions Y are accurately reproduced. 

The free pulse response of linear, time-invariant, dis- 
crete systems is given by a function known as the 
Markov parameter, 

Y m [n] = CA n ~ 1 B . (26) 


The superposition principle states that a system re- 
sponse to any arbitrary input can be obtained from 
a linear combination of impulse responses from that 
system. The generalized Hankel matrix of impulse 
responses is related to the Markov parameter by the 
superposition principle. The Hankel matrix is formed 
by windowing the impulse response data. A total 
of K data points are provided at discrete time steps 
n = and the r x s matrix H rs is formed as 

follows, 


Y m [n] 
Y m [1 + n\ 


Y m [n + s - 1] 

Ym [1 + U + S — 1] 


(27) 


1 Ym [r - 1 + n] 


Y m [r - 1 + n + s - 1], 


where s is the total size of the data window, and r is 
the number of time steps used to shift the data win- 
dow. The choice of r and s is arbitrary as long as 
r + s + n<K- 1-2. 

The ERA method eliminates redundant data by us- 
ing Singular Value Decomposition (SVD) on P° s , 


where 0 P and 0 m as the null matrices of order p and m 
respectively, and I p and I m are the identity matrices 
of order p and m. 

Since the discrete time step At = tk+i — tfc is con- 
stant, the continuous form of the discrete state-space 
realization (Eqn. (25)) is easily obtained. The contin- 
uous from, shown below, may require additional state 
dimensionality when the discrete realization has real, 
negative poles: 

X(t) = AX(t) + Bu(t ) 

Y(t) = CX(t) . (31) 

Aeroelastic ROM Development 

The full-order vector of fluid variables U(t) repre- 
sents the spatially discretized flow field obtained from 
the full-system flow solver. POD provides a trans- 
formation <I> that maps U(t) to a low order vector of 
modal coefficients U(t) (from Eqn. 18a). The reduced- 
order fluid variable U(t) will be denoted Y f, which is 
the vector of outputs Y{t) (Eqn. (31)) . 

A state-space model for Y f can be obtained from 
impulse responses using the ERA method. Impulses 
for the fluid system use the plate position and ve- 
locity coefficients (Y s ) as the forcing terms. Each 
structural term is impulsed, and the fluid system re- 
sponse is generated using the full-order model. The 
time history of the impulse response is projected onto 
each of the POD basis functions to obtain the impulse 
response of the reduced-order fluid vector Y f. POD 
basis functions are obtained using the method of snap- 
shots as described previously. The process is repeated 
for each structural mode, and the collection of impulse 
responses is used to generate a linear state-space model 
for the reduced order fluid system, 


H° rs = PDQ t . (28) 

Unwanted state dimensionality is eliminated by trun- 
cating the elements of P, D, and Q associated with 
very small singular values of H® s . The number of 
states is reduced to a minimal number q. The number 
of observations p and the number of forcing terms m 
are known from the problem formulation. The dimen- 
sion of the Markov parameter Y m [?r] is p x m. Algebra 
is used to recast Eqn. (26) in terms of the time shifted 
Hankel matrix H} s , and the elements P, D 1 and Q. 
The state-space realization flows from this manipula- 
tion, and is as follows: 


A = 

D- 1 iP T H 1 rs QD- l i , 

(29a) 

B = 

D^Q T E m , 

(29b) 

C = 

PjPP»3 . 

(29c) 

Ep and E ^ are defined below: 



— • • • j Op] , 

(30a) 

El = 

\^mi 0m 5 • • • ? 0m] 5 

(30b) 


X f = A f Xf + B f Y s , (32a) 

Y f = C f X f , (32b) 

where X f is the state vector from the ROM realiza- 
tion, and Y s represents the modal coefficients for the 
structural deformation. 

The structural model from Eqn. (7b) is coupled to 
the reduced-order fluid model (Eqns. (32a) and (32b)) 
through the projected pressure term P . The reduced- 
order variables, Y s , are obtained from the dynamic 
states X f by the linear mapping C / . Fluid pressure on 
the panel is extracted from Y s using the portion of the 
reduced order mapping (Eqn. (18a)) that pertains to 
the moving boundary. Equation (5) is used to project 
the pressure values into P. The mapping of reduced 
order fluid states to projected pressure is denoted fp, 

P = f P (X f ) . (33) 

Equation (33) is used to couple the structure and fluid 
dynamic state variables, 

Y s =R s (Y a ,fp(X f );n,X,h/L) . (34) 
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Equations (32a) and (34) comprise the low order, 
aeroelastic ROM with linear fluid dynamics, and non- 
linear plate dynamics. 

The fluid and structural terms are grouped into ar- 
rays Y and R as follows: 


Y 

R 


Xf 

Y ’ 

1 s 


A f X f + B f Y s 
R S (Y s , f P (X f ); /z, A, h/L) 


(35) 

(36) 


The reduced-order, aeroelastic system is denoted as 
simply, 


Y = R(Y) . (37) 


this fashion since any forcing function can be assem- 
bled from a series of impulses. 

The supersonic flow field was essentially linear 10 
once LCO was fully developed. Shock waves formed 
at the ends of the panel, and although they varied in 
strength dynamically with the LCO, they were station- 
ary, and the flow field between the shocks (directly over 
the panel) was linear for this case. However, impuls- 
ing the uniform, supersonic flow field with structural 
modes (and modal velocities) produced very nonlinear 
transient behavior. Figure 1 shows density contours of 
the flow 1.10802 nondimensional time units from the 
impulse of the fundamental structural velocity mode. 
The sudden appearance of a velocity profile on the 


Time Integration of the Aeroelastic System 

The aeroelastic ROM (Eqn. (37)) is integrated in 
time with the two-time level, second-order accurate, 
Crank-Nicolson method: 


yn+l _ y-n 

A t 


(i ? n+1 


R r 


(38) 


1Z = F” +1 - ~ yn ~ = °- (39) 


At each time level, y = Y n+1 is computed from Eqn. 
(39) using a chord technique with a time- frozen Jaco- 
bian 


( 7 ' T io ) (■ yk+1 ^ yk ) = -K{y k ) , (40) 

where k is a subiteration index and J a is the Jacobian 
of the reduced order aeroelastic system, evaluated for 
the base flow condition and Y s = 0. A suitable num- 
ber of subiterations are computed at each time step to 
obtain a good approximation to Y n+ ; typically, 1-2 
subiterations are generally sufficient to drive 7 Z to near 
machine zero. Since peak panel deflection is no more 
than 2% of panel length for the cases considered, the 
chord method is rapidly convergent. Prior to subiter- 
ation, y is predicted from the explicit formula 

y = Y n + AtR n . (41) 

Results 

The results that follow consider supersonic free 
stream flow conditions at Mach 1.2, with sea level con- 
ditions. The Galerkin panel model contains 4 modes, 
for a total of 8 DOFs. 

Impulse Response of a Supersonic Fluid 

Impulsing the forcing term (or terms) of a truly lin- 
ear system produces a response that is the building 
block necessary to recreate the system output from 
any arbitrary forcing function. Linear superposition 
allows the response of the system to be constructed in 


2 

1.5 

> 1 

0.5 

0 

2 

1.5 

> 1 

0.5 


Fig. 1 Density contours 1.10802 time units after 
impulse 

boundary (and its sudden removal one time step later) 
produced a shock wave of varying strength running the 
length of the panel. Early during the transient period, 
this shock welled upward away from the panel, and 
converted downstream. The convection of the shock, 
combined with the patterns of varied intensity pro- 
duce the odd (but physical) spatial oscillation above 
the panel shown in Fig. 1. After about 2.5 time units, 
this pattern had both converted well beyond the end 
of the panel, and diffused into a more benign flow pat- 
tern. After 25 time units uniform flow was restored. 
The flow dynamics were essentially linear after the ini- 
tial 2.5 time unit transient. 

The results that follow will detail the usefulness of 
such impulse responses for generating a reduced-order 
fluid model. The linear portion of the impulse response 
time integrations contained some of the LCO flow field 
characteristics, otherwise none of the ROM options 
would have reproduced LCO. However, the impulse 
responses themselves would not produce POD basis 
functions capable of correctly modelling the LCO flow 
field. This suggests that the supersonic flow field was 
not truly linear. 



Level 

Density 

10 

1.00011 

9 

1 .00009 

8 

1 .00006 

7 

1 .00004 

6 

1 .00002 

5 

0.999997 

4 

0.999975 

3 

0.999953 

2 

0.99993 

1 

0.999908 


-- 'a ’ 


k /-Vrf 
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Identification of Fluid Modes 

Fluid modes were obtained using POD as outlined 
previously. Aeroelastic fluid modes were obtained 
from a set of 100 snapshots. Snapshots of the full- 
order, aeroelastic system were taken at equally spaced 
intervals, from start-up through 25 time units, using 
time integration of the full system, with A = 25. For all 
time integration cases that follow, the aeroelastic sys- 
tem was initialized with the base flow condition, and 
a small perturbation (height of 0.0001) in the funda- 
mental panel position mode (denoted in Eqn. (7a). 
Snapshots were taken of primitive fluid variables p, u, 
v , and Et , not the conserved variables given in Eqn. 
(8b). Primitive variables enable Galerkin projection 
for the fluid as a possible means for obtaining non- 
linear terms 13 in future analysis. At Mach 1.2, LCO 
required about 300 time units to become fully devel- 
oped, and the small, 25 time unit training window was 
shown to be adequate in previous work. 10-12 

The results that follow refer to two cases. For the 
first case, the base flow term Uq from Eqn. (18a) con- 
sisted of uniform, free-stream conditions everywhere 
throughout the domain (referred to as slug flow). The 
second case considered steady state flow over the ini- 
tial panel perturbation as the base flow term Uq. For 
both cases, the first two modes for each fluid variable 
contained over 98% of the energy content, and sys- 
tem realization was performed using a total of M = 8 
fluid modes (2 modes per fluid variable). We also at- 
tempted to use the impulse response data as snapshots, 
in lieu of aeroelastic time integration. The same sys- 
tems realization procedure was repeated to produce 
a Volterra-POD ROM for this third case, but this 
ROM did not correctly produce LCO. We documented 
our observations, and recommendations regarding this 
third approach in a separate section. 

System Realization 

We considered 8 POD modes, the dimensionality of 
Y f, which produced eight impulse responses for each 
forcing function. With 8 forcing terms in Y s , the total 
number of impulse responses numbered 64. Realiza- 
tion via the ERA process for each of the aeroelastic 
cases is detailed below. 


data was windowed using s = 192 and r = 100. Every 
fifth data point was used in the realization algorithm 
providing data at the rate of dt = 0.07716. These val- 
ues of dt, s and r were chosen by trial and error to 
produce realizations whose impulse responses closely 
matched the data. The value of m was 8 to match 
the number of forcing terms, and p = 8 was chosen to 
match the number of ROM coefficients. The collection 
of impulse responses formed an 8 x 8 Markov parameter 
Y m [n] function from Eqn. (26). The number of states, 
q — 8, was chosen to match the number of ROM co- 
efficients, so SVD on the Hankel matrix formed from 
Y [n = 1] was used to truncate all but the largest 8 
singular values, yielding the matrices P, D , and Q. 
Equations (29a, 29b, 29c) were then used to generate 
a linear state-space model for the reduced-order fluid 
system, 

x[n+l\ = Ax[n\ + B (3 u[n] , (42a) 

y[n] = Cx[n ] , (42b) 

where (3 was a scaling parameter that was used to cal- 
ibrate the forcing amplitude. 

The value of the scaling parameter (3 = 650 was 
set by tuning the ROM results to the snapshot data. 
Theoretically, (3 should have been the inverse of the 
impulse size {(3 = 1/0.1 = 10). The need for an or- 
der magnitude increase in (3 reveals an inefficiency in 
projecting the impulsed flow-field onto the aeroelastic 
modes. Evidently, the impulsed flow-field contained 
structures not adequately represented in the aeroe- 
lastic modes, and a significant amount of flow energy 
was not captured in the projections used to compute 
the modal-impulse behavior. However, enough linear, 
aeroelastic information was resident in the impulsed 
flow-field for the Volterra-POD realization to produce 
correct results (with (3 properly adjusted). 



Uniform Base Flow 

A state-space realization of the form in Eqn. (25) 
was obtained using ERA for the slug-flow base case. 
The impulse amplitude was arbitrarily chosen to be 
0.1. The full-system response to this impulse was sam- 
pled over 30 non-dimensional time units at a rate of 
dt = 0.015432 for a total of I< = 1944 discrete data 
points. The fluid system impulse response was gen- 
erated using the full-order model. The time history 
of the full-order impulse response was projected onto 
each of the POD basis functions to obtain the impulse 
response of the reduced-order fluid variable Y f. The 



Fig. 2 Response of density modes to velocity term 
impulse, uniform base flow case 

The impulse response of Eqns. (42a and 42b) was 
obtained with (3 = 1. The impulse responses of the 
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reduced-order system were in good agreement with the 
impulse responses from the full-order system used by 
ERA. The response of the two density modes to an 
impulse in the first panel velocity term b \ within Y s 
(defined in Eqn. (7a)) are shown in Fig. 2. The con- 
tinuous form of Eqns. (42a and 42b) were obtained 
via a function call in MATLAB, which provided the 
matrices Af, Bf and Cf for time integration of the 
aeroelastic model given in Eqn. (37). 

Steady-state Base Flow 

The same procedure was repeated, but the aeroe- 
lastic modes were computed using a steady-state base 
flow condition described earlier. For this realization, 
every 20“ data point was used in the realization algo- 
rithm providing data at the rate of dt = 0.3086. The 
same impulse amplitude of 0.1 was used for this case, 
and the full-system impulse response was sampled at 
the same rate. However, the impulses were added to 
the steady-state panel deflection for this case. The 
data was windowed using s = 47 and r = 20. Again, 
the value of m was 8 to match the number of forcing 
terms, and p = 8 was chosen to match the number of 
ROM coefficients. The scaling parameter from Eqn. 
(42a) was (3 = 800, and was determined by tuning the 
ROM to the snapshot data. 


ROM Time History 

Both the slug-flow base case and the steady-state 
base case Volterra-POD ROMs, described above, were 
time-integrated using the aeroelastic training condi- 
tions of Mach 1.2 and A = 25. The results are shown 
in Fig. 4. Both cases correctly predicted LCO, but 




Fig. 4 Panel deflection ( Wd/h ) time history, A = 25, 
Mach 1.2 



Fig. 3 Response of density modes to velocity term 
impulse, steady-state base flow case 


The larger value of dt eliminated much of the high- 
frequency transient, and focused ERA on the low- 
frequency portion of the impulse response. This is 
reflected in the impulse response accuracy shown in 
Fig. 3, which considers the response of the two density 
modes to an impulse in the first panel velocity term b\ 
within Y s (defined in Eqn. (7a)). The first 3 seconds 
of the impulse response curve was not matched very 
closely by the ROM; however, these initial transients 
were not important to the LCO flow field. 


the steady-state base case was more accurate in am- 
plitude, frequency and phase than the slug-flow base 
case. The 3/4 chord panel-amplitude was muted by 
9% for the steady-state base case, and magnified by 
25% for the slug- flow base case. The phase and fre- 
quency error were negligible for the steady-state base 
case, but the larger amplitudes of the slug-flow base 
case introduced a small increase in LCO frequency (re- 
sulting in an accumulating phase error). 

We suspect the improvement in performance asso- 
ciated with the steady-state base flow was most likely 
due to the choice of data windowing parameters used 
in the ERA realization. We noticed no substantive 
differences in either the POD modes, or the impulse 
response of the full system flow-field between either 
case. Data windowing parameters were selected to 
provide a realization whose impulse response closely 
matched the original impulse data. As noted previ- 
ously, reducing the size of dt introduced high-frequency 
dynamics into the realization. The slug-flow base case 
was formed using a very small value of dt. Conse- 
quently the impulse response of the model (see Fig. 2) 
matched the initial transient in the data better than 
the steady-state base flow case (i.e. Fig. 3), which 
used a much larger dt. However, the high-frequency 
data in the impulse response was not germane to the 
large-time behavior of the aeroelastic system, and its 
inclusion resulted in a less accurate ROM under aeroe- 
lastic conditions. 
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ROM Robustness 

Both aeroelastic ROMs were time-integrated using 
a variety of panel dynamic pressure values (A). The 
intent was to explore the predictive accuracy of the 
Volterra-POD ROMs across a nonlinear parameter 
space. Both ROMs were trained at A = 25 (as de- 
scribed earlier), and robustness was defined as the 
ROM’s ability to predict panel amplitude (at the 3/4 
chord position) in fully developed LCO across the pa- 
rameter space, including non-LCO cases. 



Nondimensional Dynamic Pressure, X 


Fig. 5 Panel response verse dynamic pressure 

The ROM results are compared with full system re- 
sults, and results from the literature 30 in Fig. 5. The 
steady-state base case ROM was better suited to large 
LCO amplitudes where the larger panel amplitudes ex- 
cited a panel nonlinearity that corrupted the results 
of the slug- flow base case. Conversely, the use of slug 
flow as the base-flow term permited a more accurate 
prediction of LCO onset at the lower values of A. A 
non-LCO solution for the slug-flow base case occurred 
when Y f = [0] and Y s = [0] , but the use of steady- 
state flow over an initial, non-zero value of Y s required 
a non-zero value of Y f to produce Y a = [0]. The 
steady-state base case had difficulty producing this 
result. In addition, the small panel amplitudes near 
LCO-onset did not excite the high-frequency errors in 
the slug-flow base case that were evident at larger val- 
ues of A. 

Aeroelastic Structural Modes 

Aeroelastic structural modes were generated using 
100 snapshots of the structural response obtained dur- 
ing the training of the fluid ROM. The structural 
snapshots corresponded exactly in time with the set 
of 100 snapshots used to construct the fluid ROMs. 
Snapshots were taken of the panel position and ve- 
locity vectors (w(t) and w(t) respectively), and sub- 
space projection 8 was used to form a reduced-order 


structural model. Subspace projection relied on the 
Galerkin panel model for time integration. The panel 
position and velocity from the Galerkin panel model 
were projected onto the POD basis functions at ev- 
ery step in the time integration. Subspace projection 
demonstrated the adequacy of the POD modes at cap- 
turing the dynamic panel behavior, while maintaining 
the nonlinearity of the panel dynamics. 




Fig. 6 Panel response ( Wd/h ) with aeroelastic 
structural modes 

The reduced-order structural model was tightly cou- 
pled with the steady-state base case fluid ROM and 
time-integrated using the parameter value A = 25. 
Figure 6 compares the results with the full system 
response, and the POD/ROM results from Fig. 4 
(for the steady-state base case). For clarity, the top, 
right-hand corner of the entire time response is ex- 
panded in the lower portion of Fig. 6. Two POD 
modes per structural variable (4 DOFs total) yielded 
essentially identical results to the 4 mode (8 DOF) 
Galerkin result. Further order-reduction greatly de- 
creased the panel response. For this problem, the 
full-system structural model was very low order, and 
the additional order reduction from POD was immate- 
rial. However, future application of this technique will 
involve very high-order, nonlinear structural models 
requiring order-reduction along with the fluid model. 
These results demonstrate that a single training event 
can produce adequate POD modes for both the fluid 
and structure. 

ROMs Using Impulse Response Modes 

As an excursion, the aeroelastic fluid modes were 
replaced with POD modes derived from the impulse 
responses of the full system. Forty snapshots were col- 
lected from each of the eight impulse responses gener- 
ated by impulsing the elements of Y s . The snapshots 
were taken at even intervals over a time integration 
lasting 30 time units. Since there were 8 impulse re- 
sponses, a total of 320 snapshots were generated. Over 
98% of the flow energy was contained in the first 4 
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modes for each of the 4 fluid variable. This suggests a 
16 mode ROM would be adequate for the fluid. To test 
the adequacy of the POD basis, the full-system flow 
solution from time integration using A = 25 and Mach 
1.2 was projected onto the POD basis during time in- 
tegration. The modal amplitudes were examined and 
modes 6, 7 and 8 in energy content had much greater 
contribution to the aeroelastic flow field than did the 
first 4 modes (note that the modes were normalized). 
This observation led us to consider 10 modes per fluid 
variable for ROM realization. ERA was used to syn- 
thesize a 40 mode ROM, which required 43 states after 
conversion to the continuous form. The Volterra-POD 
ROM was coupled with the Galerkin panel model as 
described previously. 



Time 


Fig. 7 Panel pressure degradation using impulse 
response modes 

Unfortunately, the aeroelastic ROM produced a 
slowly divergent flutter instead of LCO. However, the 
panel pressures from the Volterra-POD ROM corre- 
sponded very well with the pressures obtained by pro- 
jecting the full-system response onto the POD basis, 
and extracting the panel pressures (see Fig. 7). This 
further suggests that the POD basis derived from the 
impulse response data was not adequate for modeling 
the aeroelastic flow field. The accuracy of the ROM 
at representing the full system response was encour- 
aging, and the frequency of unstable oscillations did 
match the LCO frequency. There is clearly some flow- 
field structure in the impulse response modes relevant 
to LCO, but more investigation is required before these 
modes can be used to generate an accurate ROM for 
the supersonic case. 

Computational Performance 

The motivation for employing the POD-Volterra ap- 
proach was to realize a computational performance 
improvement consistent with the reduction in the num- 
ber of DOFs. Computational performance, summa- 


rized in Table 1, was assessed by measuring the wall 
clock time to provide fully developed LCO of the flow 
(400 time units with A = 25). All computations were 
run on a 800 MHz Pentium-based, personal computer. 
The time-step used for the full system was based on 
a Courant-Friedrichs-Lewy (CFL) condition of 0.5, 
which was the highest value allowing for stability of 
the second-order method. The reduced-order time in- 
tegration used a time step size of 0.05 time units. The 
POD-Volterra ROM reduced compute time by four- 
orclers-of-magnitude, and realized an improvement in 
performance consistent with the DOF reduction. 


Flow Solver 

Fluid DOFs 

Wall Clock Time 

Full-order 

64400 

22304 sec 

POD-Volterra 

8 

8.132 sec 


Table 1 Computational performance 


The cost of computing the state-space realization 
using ERA was small relative to the cost of a full- 
system analysis. Since there were eight forcing terms, 
eight additional runs of the full order solver were re- 
quired to provide the impulse response data. Each 25 
time unit impulse response runs required 1394 seconds 
wall-clock time, resulting in 11152 seconds of computer 
time for generating the impulse response data. An ad- 
ditional 25 time unit run was required for snapshot 
collection to generate the aeroelastic modes, bringing 
the total computer processing time to 12546 seconds, 
or roughly half of the computational cost associated 
with a single full-system time integration. 

Conclusions 

The Volterra-POD approach produced a stable and 
accurate aeroelastic ROM with four orders of magni- 
tude reduction in problem order and computational 
expense. Aeroelastic modes were formed from snap- 
shots obtained during the initial build-up of LCO, and 
with a fixed value of dynamic panel pressure. The 
full system model was impulsed using the structural 
velocity and position modes. The flow field was pro- 
jected onto the aeroelastic modes to determine the 
modal amplitudes, and a linear, state-space realiza- 
tion for the fluid dynamics was synthesized from the 
modal impulse responses. The fluid and structural 
models were tightly coupled to form the aeroelastic 
ROM. Two cases were considered, one used uniform 
flow as a base term about which perturbations were 
computed by the reduced-order fluid model. The sec- 
ond case used steady-state flow over the initial panel 
deflection as the base-flow term. Both cases resulted 
in ROMs that correctly predicted LCO behavior over 
a wide parameter space; however, we concluded that 
uniform base flow was more desirable for predicting 
LCO onset. In addition, it was desirable to filter high- 
frequency information out of the impulse response data 
by computing the realization with relatively large time 
steps. 
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Fluid modes obtained from snapshots of the im- 
pulsed flow-field were found to be inadequate for accu- 
rate reduced-order modeling of the supersonic, aeroe- 
lastic flow field. However, we only made an initial 
attempt to examine this approach, and additional re- 
search should be conducted before dismissing the tech- 
nique completely. 

Reduced-order modeling of the structure was also 
explored, and the results demonstrated that a single 
training event could produce adequate POD modes for 
both the fluid and structure. Future applications of 
this technique will involve very high-order, nonlinear 
structural models requiring order-reduction along with 
the fluid model. 
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